A new algorithm for numerical simulation of Langevin equations 
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Formulated is a new systematic method for obtaining higher order corrections in numerical simulation of 
stochastic differential equations (SDEs), i.e., Langevin equations. Random walk step algorithms within a given 
order of finite At, are obtained so as to reproduce within that order a corresponding transition density of the 
Fokker-Planck equations, in the weak Taylor approximation scheme E). A great advantage of our method is its 
straightforwardness such that direct perturbative calculations produce the algorithm as an end result, so that the 
procedure is tractable by computer. Examples in general form for curved space cases as well as flat space cases 
are given in some order of approximations. Simulations are performed for specific examples of U(l) system and 
SU(2) systems, respectively. 



1. Introduction 

We propose a new systematic algorithm for ob- 
taining higher order corrections in numerical sim- 
ulation of SDEs, i.e., Langevin equations in either 
flat or curved space. 

The problem is how to realize random walk step 
algorithms for At, so as to reproduce proper cor- 
rections within a given order 0(At n ) of corre- 
sponding transition density which evolves accord- 
ing to the Fokker-Planck equation. 

Our systematic algorithm consists of all di- 
rect perturbative calculations within some order 
0(At n ) such as normal ordering of the time evo- 
lution operator, exponentiation of polynomials, 
completing squares of polynomials, and inversion 
of polynomials. These direct perturbative calcu- 
lations produce the random walk step algorithm 
as an end result without any matching equations, 
so that the procedure is tractable by computer. 

It is to be noted that Drummond et al.0] as- 
sumed linear step algorithms with respect to ran- 
dom variates, but we do not. In their scheme, 
equations becomes overconstrained and will have 
no solutuion in higher orders, since otherwise the 
obtained algorithm generates Gaussian distribu- 
tion for the transition density, which is, however, 
generally impossible in higher orders. 



2. Formalism 

The Langevin equation in d dimensional curved 
space with an initial condition x — xq at t = to, 
reads as 

fa* =u l (x)+e* a (x)ri a (t) (i, a = 1, 2, • • • , d),(l) 
<77 Q (^(t') >=2a al3 5(t~t'), (2) 
and 

<V a (t)>=0, (3) 

where (1) is understood as Ito's SDE. Then we 
can derive the corresponding Fokker-Planck equa- 
tion for the transition density in the curved space 
with an initial condition P(to, x; to, xq) = S(x — 
xq) as 

d t P(t, x; to, x Q ) = KP(t, x; t , x ), (4) 
where 

K = d i (d j g i i(x)-u i (x)), (5) 
and 

g ij = eW"' 3 - (6) 

Covariant form of the Fokker-Planck equation 
can be derived by identification, cj) = PJ -J detgij 
and in case that 

u\x) = (d jg »(x)P (x))/P (x), (7) 
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P(x) tends to a unique asymptotic distribution 
Po(x). 

Given a u l {x), we are to find a random walk 
step algorithm Ax* , a function of random variates 
r/ a , so as to satisfy within a given order 0(At n ), 

P(t + At, x- t , x ) =< 5{x - (x + Ax)) > (8) 

where 

< rfrf >= 25 a(3 At, (9) 
and 



where z l is a polynomial of z\z % — z % {\ft, z) such 
that z 1 — z % + c 1 + ■ ■ ■, and TV is a normalization 
constant independent of z. Using the 6— function, 
we can write (15) as 



P(t + At,x;t ,x ) 



N 



V(47rAt) d 



x / drje~^ 6 d {z i -6 i ), 



(16) 



< T] a >= 0. 



(10) 



where l = \/g^{xo) ff ' , and inverting the argu- 
ment of the <5— function with respect to z l , we 
rewrite the 5— function as 



2.1. A new systematic algorithm for higher 
order random walk step algorithms 

First we write a truncated exponential operator 
of the 0(At n ) formal solution 



s AtK S{x - x ) 



5 d {~z l - 8 l ) = det{dz/d6)5 d {z l - z l {6)) 
to obtain 

N 



(17) 



N f 

(11) P(t + At, x; t ,x ) = -. xj / dfjdet(dz/d6) 

\ 4tt ) d J 



in normal order, i.e., differential operators in the 
leftmost position, and then all velocity and met- 
ric fields and their derivatives may be evaluated 
at the position x and become commutable with 
differential operators. Then we insert 

e -AtgV {x )d i d j e Atg^ (x )d i d j 



xe-i f ' 2 S d (x i -xi-Ax^O)), 



(18) 



where Ax 1 = \/Atz t {9). The final part of our 
procedure to reach the goal in the form 



between e AtK and S, and make use of 



5 Ats ' 3{xo)9 ^6(x-x ) 



^/(4irAt) d g 



(12) 



P(t + At,x;t ,x ) 



1 



dr/e 4V 



e 4 



where z l = (x l — x l )/V At, g = det(g %J (x )) , and 
di = d z z /\/Ai to obtain the transition density in 
a form as 



,(13) 



TP 

xd d (x l - x - Ax l (8)), (19) 
is to find a transformation ff — ff {rf) such that 

dfjdet(dz/d6)e-i i > 2 = 1 drje' (20) 



P(t + At, x; to,x ) = F(v At, z, ■ 



V(4^Ai)% 



e 4i 



(14) 



which we prove to be always possible as follows. 

Algorithm to find the transformation from ff 
to rf: Noting that 



det{dz/d6) = ! + ■■ 



(21) 



We use a singular coordinate z % = [x % — 
x l )/y/At due to a singular nature of the pertur- 
bation expansion. 

Within the given order, we exponentiate the 
polynomial factor F(\* r Ai,z,- ■ ■), and proceed 
perturbatively order by order so as to elimi- 
nate the leading order corrections in the expo- 
nent polynomial by finding proper transforma- 
tion, that is, by completing squares, to obtain 

P(t +At,x;t ,x ) = n = = e-^*°>^,(15) 
v /(47rAt)% 



we exponentiate the determinant factor, and con- 
sider its exponent in a form of polynomial of fj as 
— l/4(ry 2 + •••). We assume that constant terms 
are absent since they can be factored out as nor- 
malization constant, and also that either linear 
or quadratic terms with respect to fj in the low- 
est leading order correction terms are absent since 
they can be eliminated by some linear transforma- 
tion which only effects a constant Jacobian fac- 
tor. Let us suppose that the least power terms 
in the leading order corrections may be of a nth 
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power of fj, then we can eliminate it by a suitable 
transformation, i.e., by completing square, which 
simply produces a (n-2)th power term in the ex- 
ponent when the corresponding Jacobian factor is 
exponentiated. This shows that the leading order 
terms can be transformed away by finite steps to 
give the r.h.s. of (20), and completes our system- 
atic approach. 

For the sake of practical validity of the method, 
we applied, in use of suitable computer softwares, 
the above algorithm to cases in the flat space and 
in the curved space within some orders 0(At n ), 
respectively. The results are new as far as we 
know. 

2.2. Flat space local order 3 algorithm 

In flat spaces i.e. g^ v = <5 M „ the Langevin step 
algorithms that reproduce the transition density 
of the Fokker-Planck equations up to order At 3 
are obtained in the Taylor scheme and in the 
Runge-Kutta-like scheme, respectively. In the 
latter the higher order derivatives are absent. 

Our final result in the Taylor scheme is 

Ax\s, 77, x ) = srf + s 2 u l + ^s 3 rfdjU l 

+s i (lv 3 V k d J d k u i + -uidjv? + ~c>V) 

+s\\rfu k d i d k u i + Irfd^dkU* 
3 6 

+s\\d 1 u k d 1 d k u l + lv?u k d i d k u i + -5V0,-u* 
3 6 6 

+ \u 3 d 3 u k d k u l + -u'fyflV + \d 2 d 2 u l ) (22) 

D o 

Here s = \/At and 77* 's are the Gaussian random 
variates which satisfy 

//'.//' LM"< (23) 

Remark that there are terms non-linear in 77's. 

We also derived the 0(At 3 ) Runge-Kutta-like 
[] scheme which contains two sets of the Gaussian 
random variates rf and r}'° . 

1 Transcription of (22) to the pure Runge-Kutta scheme 
is impossible due to contraction structure of rf and u J in 
(22). 



2.3. Curved space local order 2 algorithm 

In general positive metric curved spaces, the 
Langevin step algorithm that reproduces the 
transition density of the Fokker-Planck equation 
up to order At 2 is obtained in the Taylor scheme. 

Our final result is 

Ax l (s 1 77, xq) = s9 % 

+ s 3 ( + U k ^d jg lk + - A e kg mn d m d n9 lk 

~0j9 km 9nid k g in d m g jl - ^jir^/) 
+s\^d J u t + ^ k djd k u l ) (24) 

Here all g lJ 's, u l, s and their higher derivatives are 
evaluated at xq, and 9 z, s are given by 

0* = VgWo)v 3 (25) 

and 77j's are the Gaussian random variates as in 
(23). 

As an example, we transcribe the above al- 
gorithm to the standard SU(2) algorithm on S 3 
manifold. The S 3 manifold can be patched by two 
hemispheres and points x on the northern hemi- 
sphere are projected stereographically at y on the 
tangential plane at the north pole and points on 
the southern hemisphere are projected on the tan- 
gential plane at the south pole. 

The ^-component of a point on the northern 
hemisphere is parametrized as X4 = cos8 and y = 
2ton||^y . Similar parametrization is done for the 
southern hemishere. 

Using 3 dimensional Gaussian random variates 
rf of variance 1, the Ay 1 can be written as 

Ay^s, 77,2/0) = aVUfif 
+s 2 [-(Af3 + f)f + ^f(y k r 1 k )/2] 

+s 3 [rf{- 12/3/ + 8(3 +I(/ 2 + /)) 

-f/(i/vy]/>/2 

+s iV -m 2 + 4(2/ - 1)(3 -\{f + /)] (26) 
where / = 1 + y 2 /4. 
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The expectation value of X4 is 



SU2 S 3 [3=1 



< X4 > 



1 N 



AT * — ' y z 

i=\ J 



(27) 



3. Numerical examples 

We consider the standard probability distribu- 
tion 

p(x) cx exp(/3S(x)) (28) 

where x — 9 and 5* = cosO in U(l) case, and 
S = 4x4 with a four-vector Xi satisfying xj = 1 
in SU(2) case, respectively. 

3.1. U(l) on S 1 

The numerical simulation of the expectation 
values < sin 2 9 > is performed, and its data are 
shown in Fig.l. The data are taken from succes- 
sive 2000 Langevin steps each of which are the 
average of 100,000 runs. 
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Figure 1. The results of simulation of < sin 2 9 > 
in U(l) system. (3 = 5, standard method in the 
Taylor scheme. 




0.05 0.1 0.15 0.2 
At 



Figure 2. The results of simulation of < x 4 > of 
SU(2) in S 3 . [3 = 1 



4. Conclusion and Outlook 

General algorithm for obtaining higher order 
corrections to the Langevin step algorithms with 
respect to At was formulated. The method is 
straightforward without matching equations and 
is tractable by computer. 

A local third order algorithm for simulating 
general flat space Langevin equations and a lo- 
cal second order algorithm for simulating general 
curved space Langevin equations were obtained. 

This method is based on Gaussian distribution 
in use of noncompact coordinates, as its lowest 
leading order approximation to transition den- 
sity. Even in case of compact curved spaces as 
S™, however, patching algorithm in use of stereo- 
graphic projection on two noncompact tangential 
planes turned out to be successful and fast on 
computer. Similar method is applicable in prin- 
ciple to simulation of imaginary time evolution of 
wave functions in quantum mechanics. 
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